%% fig 2a  population vector correlation matrix example 

savepath='/media/Big1/15tracks/placefieldCorr/';

load(fullfile(savepath,'PlFieldsCCMatrixAlldir_geometryFeatures.mat'),'res','datacc',...
    'm','se','pval');

ires=2
trackx=res(ires).trackx1;
plotind=trackx(:,1)<=15 & trackx(:,2)==1;
ccr=res(ires).ccr2dirmax(plotind,plotind);

xs=[];ys=[];
for i=1:size(ccr,1)
    for j=1:size(ccr,2)
   if i>=j;ccr(j,i)=NaN; xs=[xs;i];ys=[ys;j];end
    end
end

f1=figure;%colormap copper
imagesc(ccr);colorbar('Ticks',0:.2:1);
hold on;
plot_square_on_imagesc2(xs,ys,'facecolor','w','edgecolor','w');

xlabel('Track');ylabel('Track');title('Correlation');
set(gca,'xtick',5:5:15,'ytick',5:5:15);
setaxisformal

% savepdf(f1,'/media/Big1/15tracks/figures/example_pop_vec_corr.fig')
% savepdf(f1,'/media/Big1/15tracks/figures/example_pop_vec_corr.pdf')
%% fig 2b population vector correlation  --- GLM of geometry  fig1d
load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_glm.mat'),'mdl',...
    'crit_afterremove','varname','pval_eachvar','nvar');
disp(varname)
param_str={'Maze','In Out','Middle Bar','Near Curtain','Orientation',...
    'Near FixBar','Near Wall'};
plotind=1:7;
f1=figure('position',[333   134   558   367]);
bar(1:size(crit_afterremove(plotind,:),1), crit_afterremove(plotind,3));hold on;
sigstar1(1:length(plotind),pval_eachvar(plotind),crit_afterremove(plotind,3),0);
set(gca,'xtick',1:length(plotind),'xticklabel',param_str(plotind),'XTickLabelRotation',30,...
    'xlim',[0,length(plotind)+1],'box','off');
ylabel('\Delta Deviance');
setaxisformal

f2=figure('position',[333   134   558   367]);
stem(1:size(crit_afterremove(plotind,:),1), crit_afterremove(plotind,3),'o');hold on;
sigstar1(1:length(plotind),pval_eachvar(plotind),crit_afterremove(plotind,3),0);
set(gca,'xtick',1:length(plotind),'xticklabel',param_str(plotind),'XTickLabelRotation',30,...
    'xlim',[0,length(plotind)+1],'box','off');
ylabel('\Delta Deviance');
setaxisformal

%  savepdf(f1,'/media/Big1/15tracks/figures/pop_vec_geometry_glm.pdf')
% savepdf(f2,'/media/Big1/15tracks/figures/pop_vec_geometry_glm_dot.pdf')
% figure;
% subplot(221);plot_scatter([y,newf]);title(mdl.Deviance)
% subplot(222);plot_scatter([y,newf1]);title(mdl1.Deviance)
% subplot(223);plot_scatter([y,newf2]);title(mdl2.Deviance)
%% fig 2b  GLM add analogous as feature
savepath='/media/Big1/15tracks/placefieldCorr';
load(fullfile(savepath,'PlFieldsCCMatrixAlldir_glm_addAnalogous.mat'),'res','mdl',...
    'crit_afterremove','varname','pval_eachvar','nvar')
mdl

param_str={'Maze','In Out','Middle Bar','Near Curtain','Orientation',...
    'Near FixBar','Near Wall','Analogous'};

f2=figure('position',[333   134   558   367]);
stem(crit_afterremove(:,3),'o');hold on;
sigstar1(1:nvar-1,pval_eachvar,crit_afterremove(:,3)+0.01,0);
set(gca,'xtick',1:nvar-1,'xticklabel',param_str,'XTickLabelRotation',30,...
    'xlim',[0,nvar],'box','off');
ylabel('\Delta Deviance');

setaxisformal
% savepdf(f2,'/media/Big1/15tracks/figures/pop_vec_geometry_glm_dot_analogous.pdf')
%% GLM of run pairs add Direction
clear
savepath='/media/Big1/15tracks/placefieldCorr';
load(fullfile(savepath,'PlFieldsCCMatrixAlldir_glm_eachDir.mat'),...
    'mdlAll','res',...
    'dirdata','dirdata_str','nvar');
 fd=figure('position',[779   616   791   347]);
for i=1:2
    crit_afterremove=mdlAll(i).crit_afterremove;
    dirdata_str1=dirdata_str{i};
    pval_eachvar2=mdlAll(i).pval_eachvar;
subplot(1,2,i)
    stem(crit_afterremove(9,[1,2]));hold on;
    sigstar([1,2],pval_eachvar2(9));
%     set(gca,'TickLabelInterpreter','dex');
    set(gca,'xtick',1:2,'xticklabel',...
        {'Full model','Reduced model(-Direction)'},'XTickLabelRotation',20,...
        'xlim',[0,3],'box','off');
    ylabel('Deviance');
    title(sprintf('%s,p=%.3g',dirdata_str1,pval_eachvar2(9)));
end
setaxisformalAll(12)
% savepdf(fd,'/media/Big1/15tracks/figures/pop_vec_geometry_glm_eachdir_direction.pdf')
%% fig 2cef. population vector correlation  --- geometry 
clear
load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_geometryFeatures.mat'),'res','datacc',...
    'm','se','pval');

load('/media/Big1/15tracks/15trackIDnew.mat')

ntype=length(datacc);
datacc1=cell(length(datacc),length(datacc{1}));
for i=1:length(datacc)
    datacc1(i,:)=datacc{i};
end

f2=figure('color','w','position',[363    71   1380         305]);
for itype=1:ntype
    xstr={res(1).typestr{itype,1},res(1).typestr{itype,2},'diff','same'};
    subplot(1,ntype,itype)
    errorbarformal(1:4,m(itype,:),se(itype,:));
    %     hold on;
    %     for i=1:4
    %         ydata=datacc1{itype,i};
    %         plot(i+rand(size(ydata))/3,ydata,'.k');hold on;
    %     end
    sigstar({[1,2],[3,4],[1,3],[2,3]},pval(itype,:))
    set(gca,'xtick',1:4,'xticklabel',xstr,'XTickLabelRotation',30)
    ylabel('CC');
end
setaxisformalAll(10)
% 


f3=figure('color','w','position',[ 112   665   411   273]);
for itype=3%1:ntype
    xstr={res(1).typestr{itype,1},res(1).typestr{itype,2},'diff','same'};
    
    [ydata,g]=cell2group(datacc1(itype,3:4)');
    % -------- bar with dots
    %     for i=3:4
    %         ydata=datacc1{itype,i};
    %         plot(i-2+(rand(size(ydata))-0.5)/3,ydata,'.','color',0.5*[1,1,1]);hold on;
    %     end
    %
    [hdot,hbar]=plot_errorbar_with_data(ydata,g);
    
    %     plot(g+(rand(size(ydata))-0.5)/3,ydata,'.','color',0.5*[1,1,1]);hold on;
    %     he=errorbarformal(1:2,m(itype,3:4),se(itype,3:4));
    %     set(he(1),'facecolor','none','edgecolor','k','linewidth',2,'linestyle','-');
    %     set(he(2),'linewidth',2);
    % --------- boxplot
    %          boxplot(ydata,g,'symbol','','Colors','k');
    
    
    sigstar([1,2],pval(itype,2))
    set(gca,'xtick',1:2,'xticklabel',{'Orthogonal','Parallel'},'XTickLabelRotation',30)
    ylabel('CC');
end
setaxisformalAll(10)
disp(pval(3,2))

% ---- repetition of same track
sametrcc=[];
for irat=1:length(res)
    ccr=res(irat).ccr2dirmax;
    trackx1=res(irat).trackx1;
    for i1=1:size(ccr,1)
        for i2=1:size(ccr,2)
            if trackx1(i1,1)==trackx1(i2,1) && i1~=i2
                sametrcc=[sametrcc;ccr(i1,i2)];
            end
        end
    end
end
[sametrccm,sametrccse]=mean_se(sametrcc);
% ----------- between maze
betwmazecc=[];
withinmazecc=[];
for irat=1:length(res)
    ccr=res(irat).ccr2dirmax;
    trackx1=res(irat).trackx1;
    for i1=1:size(ccr,1)
        for i2=1:size(ccr,2)
            if (trackx1(i1,1)<8 && trackx1(i2,1)>7) ...
                    &&(trackx1(i2,2)==1 && trackx1(i1,2)==1)
                betwmazecc=[betwmazecc;ccr(i1,i2)];
            end
            if ((trackx1(i1,1)<8 && trackx1(i2,1)<8) || (trackx1(i1,1)>7 && trackx1(i2,1)>7)) ...
                    &&(trackx1(i2,2)==1 && trackx1(i1,2)==1) && i1~=i2
                withinmazecc=[withinmazecc;ccr(i1,i2)];
            end
        end
    end
end
[betwmazeccm,betwmazeccse]=mean_se(betwmazecc);
[withinmazeccm,withinmazeccse]=mean_se(withinmazecc);
% --------analogous
anatrcc=[];
for ii=1:length(res)
    trackidx1=trackIDx([trackIDx.rat]==res(ii).irat);
    trackx1=res(ii).trackx1;
    ia=trackx1(:,2)==1&trackx1(:,1)<16;
    trackx1=trackx1(ia,:);
    ccr=res(ii).ccr2dirmax(ia,ia);
    for i1=1:15
        for i2=1:15
            %outer track
            if trackidx1(i1).inner==0 && trackidx1(i2).inner==0
                if i2-7==i1
                    anatrcc=[anatrcc; ccr(i1,i2)];
                end
            else
                %inner track
                if i1==5 && (i2== find([trackidx1.maze]==2&[trackidx1.inner]==1&[trackidx1.nearFixBarrier]==1 ...
                        &[trackidx1.parallelCurtain]==trackidx1(i1).parallelCurtain))
                    anatrcc=[anatrcc; ccr(i1,i2)];
                elseif i1==7 && (i2== find([trackidx1.maze]==2&[trackidx1.inner]==1&[trackidx1.nearFixBarrier]==0 ...
                        &[trackidx1.parallelCurtain]==trackidx1(i1).parallelCurtain))
                    anatrcc=[anatrcc; ccr(i1,i2)];
                end
            end
        end
    end
end
[anatrccm,anatrccse]=mean_se(anatrcc);
para_notana_cc = datacc1{3,4}(~ismember(datacc1{3,4},anatrcc));
[para_notana_ccm,para_notana_ccse]=mean_se(para_notana_cc);

datacell=[datacc1(3,3:4),{sametrcc},{betwmazecc},{withinmazecc},{anatrcc}, {para_notana_cc} ];
datacellm=mean_se_cell(datacell);
[pvalout,grouppair,pval]=test_multigroup_rank(datacell);
[ydata,g]=cell2group(datacell);
xstr={'Orthogonal','Parallel','Same track','Between maze','Within maze','Analogous track','Parallel not analogous'};

f4=figure('position',[551    99   509   667]);
[hdot,hbar]=plot_errorbar_with_data(ydata,g,xstr);
sigstar(grouppair,pval);
ylabel('CC');
setaxisformal

% fa=figure;
% tmpdata={datacell{3},[datacell{1};datacell{2}],datacell{6}};
% [pvalout,grouppair,pval]=test_multigroup_rank(tmpdata);
% [ydata,g]=cell2group(tmpdata);
% xstr={'Same track','Different tracks','Analogous tracks'};
% [hdot,hbar,tmpm,tmpse]=plot_errorbar_with_data(ydata,g,xstr);
% sigstar(grouppair,pval);
% ylabel('Correlations');
% title(sprintf('p=%.3g; ',pval))
% setaxisformal


fa=figure;
tmpdata={datacell{3},[datacell{1};datacell{2}],datacell{6},datacell{2} };
[pvalout,grouppair,pval]=test_multigroup_rank(tmpdata);
[ydata,g]=cell2group(tmpdata);
xstr={'Same track','Different tracks','Analogous tracks','Parallel tracks'};
[hdot,hbar,tmpm,tmpse]=plot_errorbar_with_data(ydata,g,xstr);
sigstar(grouppair,pval);
ylabel('Correlations');
title(sprintf('p=%.3g; ',pval))
setaxisformal

fm=figure;
pval=ranksum(datacell{4},datacell{5});
[ydata,g]=cell2group(datacell([4,5]));
xstr={'Between maze','Within maze'};
[hdot,hbar]=plot_errorbar_with_data(ydata,g,xstr);
sigstar([1,2],pval);
ylabel('CC');
title(sprintf('p=%.3g',pval))
setaxisformal
% savepdf(f2,fullfile('/media/Big1/15tracks/figures/','pop_vec_geometry_eachparam.pdf'));
% savepdf(fa,fullfile('/media/Big1/15tracks/figures/','pop_vec_sametr_analogous.pdf'));
% savepdf(fm,fullfile('/media/Big1/15tracks/figures/','pop_vec_maze_samediff.pdf'));
% savepdf(f3,fullfile('/media/Big1/15tracks/figures/','pop_vec_geometry_eachparam_dot.pdf'));
% savepdf(f4,fullfile('/media/Big1/15tracks/figures/','pop_vec_geometry_eachparam_multi.pdf'));
% %% maze12 , same/diff maze
% clear
% load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_geometryFeatures.mat'),'res','datacc',...
%     'm','se','pval');
%% maze1/2 para/orth

clear
load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_geometryFeatures.mat'),...
    'res','datacc',...
    'm','se','pval');

load('/media/Big1/15tracks/15trackIDnew.mat')

ntype=length(datacc);
datacc1=cell(length(datacc),length(datacc{1}));
for i=1:length(datacc)
    datacc1(i,:)=datacc{i};
end

cc2maze=datacc1(4:5,3:4); %[maze1,maze2],[orth,para]
figure;
[data2maze,gr2maze]=cell2group(cc2maze(:));
[~,grouppair,pval]=test_multigroup_rank(data2maze,gr2maze);
plot_errorbar_with_data(data2maze,gr2maze,{'maze1-orth','maze2-orth','maze1-para','maze2-para'});
sigstar(grouppair,pval);
ylabel('place map cc');
% ------------- pairwise
load('/media/Big1/15tracks/placefieldCorr/PlFieldsCCMatrixAlldir_glm_addAnalogous.mat','res');
% x=cat(1,res.x);
% y=cat(1,res.ccr2dirmax_v);
% ana_v=cat(1,res.ana_v);
% trackpair=cat(1,res.track_pair);
% trackidmat_str=res(1).trackidmat_str;
% X=array2table(logical(x),'VariableNames',trackidmat_str);
% X.analogous=ana_v;
cc1orth=[];
cc2orth=[];
cc1para=[];
cc2para=[];

for ires=1:length(res)
    iori=find(strcmp(trackidmat_str,'parallelCurtain'));
    ind1orth=find(res(ires).x(:,iori)==0 & all(res(ires).track_pair>0 & res(ires).track_pair<8,2));
    ind2orth=find(res(ires).x(:,iori)==0 & all(res(ires).track_pair>7,2));
    ind1para=find(res(ires).x(:,iori)==1 & all(res(ires).track_pair>0 & res(ires).track_pair<8,2));
    ind2para=find(res(ires).x(:,iori)==1 & all(res(ires).track_pair>7,2));
    trackidx1=trackIDx([trackIDx.rat]==res(ires).irat);
    
    ind1orth_ana=nan(size(ind1orth));
    for j=1:length(ind1orth)
        tmp2=res(ires).track_pair(ind2orth,:);
        tmp=ind2orth( (res(ires).track_pair(ind2orth,1) ==  ...
            trackidx1([trackidx1.track]==res(ires).track_pair(ind1orth(j),1)).analogousTrack) ...
            & (res(ires).track_pair(ind2orth,2) ==  ...
            trackidx1([trackidx1.track]==res(ires).track_pair(ind1orth(j),2)).analogousTrack) );
        if ~isempty(tmp)
            ind1orth_ana(j)=tmp;
        end
    end
    ind1orth(isnan(ind1orth_ana))=[];
    ind2orth=ind1orth_ana(~isnan(ind1orth_ana));
    
    ind1para_ana=nan(size(ind1para));
    for j=1:length(ind1para)
        tmp2=res(ires).track_pair(ind2para,:);
        tmp1=[trackidx1([trackidx1.track]==res(ires).track_pair(ind1para(j),1)).analogousTrack,...
            trackidx1([trackidx1.track]==res(ires).track_pair(ind1para(j),2)).analogousTrack];
        tmp=ind2para( (res(ires).track_pair(ind2para,1) ==  ...
            trackidx1([trackidx1.track]==res(ires).track_pair(ind1para(j),1)).analogousTrack) ...
            & (res(ires).track_pair(ind2para,2) ==  ...
            trackidx1([trackidx1.track]==res(ires).track_pair(ind1para(j),2)).analogousTrack) );
        if ~isempty(tmp)
            ind1para_ana(j)=tmp;
        end
    end
    ind1para(isnan(ind1para_ana))=[];
    ind2para=ind1para_ana(~isnan(ind1para_ana));
    
    cc1orth=[cc1orth; res(ires).ccr2dirmax_v(ind1orth), res(ires).track_pair(ind1orth,:)];
    cc2orth=[cc2orth; res(ires).ccr2dirmax_v(ind2orth), res(ires).track_pair(ind2orth,:)];
    cc1para=[cc1para; res(ires).ccr2dirmax_v(ind1para), res(ires).track_pair(ind1para,:)];
    cc2para=[cc2para; res(ires).ccr2dirmax_v(ind2para), res(ires).track_pair(ind2para,:)];
end


%
figure;
xylim=[0,0.4];
subplot(121)
str1orth=num2str(cc1orth(:,[2,3]),'%d-%d');
pvalorth=signrank(cc1orth(:,1),cc2orth(:,1));
cor=copper(size(cc1orth,1)); [tmp,ia]=sort(mean(cc1orth(:,[2,3]),2));
sz=cc1orth(ia,2)*50;
scatter(cc1orth(ia,1),cc2orth(ia,1),sz,cor,'filled');
hold on; plot(xylim,xylim,'b:');set(gca,'xlim',xylim,'ylim',xylim);
text(cc1orth(:,1),cc2orth(:,1),str1orth)
xlabel('cc(maze1)');ylabel('cc(maze2)');
title(sprintf('orthogonal track pair, p=%.4g',pvalorth));


subplot(122)
str1para=num2str(cc1para(:,[2,3]),'%d-%d');
pvalpara=signrank(cc1para(:,1),cc2para(:,1));
cor=copper(size(cc1para,1));[~,ia]=sort(mean(cc1para(:,[2,3]),2));
sz=cc1para(ia,2)*50;
scatter(cc1para(ia,1),cc2para(ia,1),sz,cor,'filled');
hold on;plot(xylim,xylim,'b:'); set(gca,'xlim',xylim,'ylim',xylim);
text(cc1para(:,1),cc2para(:,1),str1para)
xlabel('cc(maze1)');ylabel('cc(maze2)');
title(sprintf('parallel trackpair, p=%.4g',pvalpara));

legend({'size:1st track','color:mean(track)'})
%% parallel and same dir/diff dir.  run pair
clear
savepath='/media/Big1/15tracks/placefieldCorr/';
load(fullfile(savepath,'PlFieldsCCMatrixAlldir_glm_eachDir.mat'),'res','mdlAll',...
    'dirdata','dirdata_str','nvar');
x=cat(1,res.x);
y=cat(1,res.ccr2dirmax_v);
trackidmat_str=res(1).trackidmat_str;
X=array2table(logical(x),'VariableNames',trackidmat_str);
idirind=strcmp(dirdata_str,'Headdir');
X.Direction=logical(dirdata{idirind});
X.plfCC=y;

orth=y(X.parallelCurtain==0);
para=y(X.parallelCurtain==1);
paradirsame=y(X.parallelCurtain==1 & X.Direction==1);
paradirdiff=y(X.parallelCurtain==1 & X.Direction==0);
   
%
f3=figure('color','w','position',[ 112   665   411   273]);
 itype=3
 xstr={'Orthogonal','Parallel','Parallel+samedir','Parallel+diffdir'};
%     [ydata,g]=cell2group([datacc1(itype,3:4),paradirsame,paradirdiff ]);
[ydata,g]=cell2group({orth,para,paradirsame,paradirdiff });
    % -------- bar with dots
    [hdot,hbar,ydatam,ydatase]=plot_errorbar_with_data(ydata,g);
[~,grouppair,pval]=test_multigroup_rank(ydata,g);
    sigstar(grouppair,pval)
    set(gca,'xtick',1:length(xstr),'xticklabel',xstr,'XTickLabelRotation',30)
    ylabel('CC');

setaxisformalAll(10)
%% fig2d. population vector correlation -- track distance 

load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_distance.mat'),'res','distance',...
    'discc','disx','discc1','x','m','se','disrcc','disccpval');
 
distance=res(1).distance;
nrat=length(res);
ndis=length(res(1).distance);
for irat=1:nrat
    disxlen=cellfun(@length,res(irat).distanceCC);
    disx=[];discc1=[];
    for id=1:ndis
        disx=[disx;res(irat).distance(id)*ones(disxlen(id),1)];
        discc1=[discc1;res(irat).distanceCC{id}];
    end
    res(irat).disx=disx;
    res(irat).discc1=discc1;
end

% ---- repetition of same track
sametrcc=[];
for irat=1:nrat
   ccr=res(irat).ccr2dirmax;
   trackx1=res(irat).trackx1;
   for i1=1:size(ccr,1)
       for i2=1:size(ccr,2)
            if trackx1(i1,1)==trackx1(i2,1) && i1~=i2
                sametrcc=[sametrcc;ccr(i1,i2)];
            end
       end
   end
end

[sametrccm,sametrccse]=mean_se(sametrcc);

% f1=figure('color','w','position',[400    74   972   446]);
% for irat=1:nrat
%     subplot(1,nrat,irat);
%     plot(res(irat).disx,res(irat).discc1,'o');
%     xlabel('Distance (cm)');ylabel('pop. vector CC');
%     title(res(irat).ratname);
%     set(gca,'xlim',[0,distance(end)+1]);
% end
% setaxisformalAll(10)

discc=cell(ndis,1);
for id=1:ndis
    for irat=1:nrat
        discc{id}=[discc{id};res(irat).distanceCC{id}];
    end
end
x=res(1).distance;
[m,se]=cellfun(@mean_se,discc);
disx=cat(1,res.disx);
discc1=cat(1,res.discc1);
[disrcc,disccpval]=corr(disx,discc1);

ranksum(sametrcc,cell2mat(discc))

f2=figure('color','w');
plot(disx,discc1,'o');
hold on;
errorbar(x,m,se,'linewidth',2);
errorbar(5,sametrccm,sametrccse,'linewidth',2);
xlabel('Distance (cm)');ylabel('pop. vector CC');
set(gca,'xlim',[0,200]);
title(sprintf('r=%.2f,p=%.2g',disrcc,disccpval));
setaxisformalAll(15,0)

savepath='/media/Big1/15tracks/figures/';



% savepdf(f2,fullfile(savepath,'trackgeometryDistanceAllrat.fig'));
% savepdf(f2,fullfile(savepath,'trackgeometryDistanceAllrat.pdf'));
% savepdf(f2,fullfile(savepath,'trackgeometryDistanceAllrat_withSameTrack.pdf'));

%% bi-directionality in repeat2
clear
savepath='/media/Big1/15tracks/HeadDirection/';
load([savepath,'head_direction_selectivity_newrand_repeat2.mat'],'dsiall')

%=========== original place field. index by place field difference
% dsi_1track=1-[dsiall.dsi_eachtrack_mean]';
% dsi_1trackall=1-cat(1,dsiall.dsi_eachtrack); dsi_1trackall=dsi_1trackall(:,1)';
%============== index by mean firing of place field
dsi_1track=1-[dsiall.dsi_eachtrack_norm_mean]';
dsi_1trackall=1-cat(1,dsiall.dsi_eachtrack_mean); dsi_1trackall=dsi_1trackall(:,1)';

bidsi=cat(1,dsiall.dsi_eachtrack);
%=================%

gooddsiind=find(~isnan(bidsi(:,1)))';

bidsi=bidsi(gooddsiind,:);

[~,grouppair,pval]=test_multigroup_rank(bidsi(:,1),bidsi(:,3));
[ccr,ccp]=corr(bidsi(:,1),bidsi(:,3));
figure;
subplot(121);
[hdot,he,m,se]=plot_errorbar_with_data(bidsi(:,1),bidsi(:,3));
sigstar(grouppair,pval);
xlabel('Repeat');ylabel('Bi-directionality');set(gca,'xlim',[0,4]);
subplot(122);
% bidsi2=bidsi;bidsi2(bidsi2(:,3)>1,3)=2;
plot(bidsi(:,3),bidsi(:,1),'o');hold on;lsline
xlabel('Repeat');ylabel('Bi-directionality');set(gca,'xlim',[0,4]);
title(sprintf('r=%.4g, p=%.4g',ccr,ccp));
setaxisformalAll

fb=figure;
datac={bidsi(bidsi(:,3)==1,1),bidsi(bidsi(:,3)==2,1)};
pval=ranksum(datac{1},datac{2});
[ydata,gr]=cell2group(datac);
[hdot,he,m,se]=plot_errorbar_with_data(ydata,gr);
sigstar([1,2],pval);
xlabel('Repeat');ylabel('Bi-directionality');set(gca,'xlim',[0,3]);
title(sprintf('p=%.4g',pval));
setaxisformalAll
% savepdf(fb,fullfile('/media/Big1/15tracks/figures/','biDirectionality_repeat2.pdf'));
%% fig 2g s2e s2f orientation selective cell 

clear 
load(['/media/Big1/15tracks/HeadDirection/head_direction_selectivity_newrand.mat'],'dsiall')
dsi_1track=1-[dsiall.dsi_eachtrack_norm_mean]';
dsi_1trackall=1-cat(1,dsiall.dsi_eachtrack_mean); dsi_1trackall=dsi_1trackall(:,1)';
%=================%

%
dsi=[dsiall.dsi]';
osi=[dsiall.osi]';
dsi_pval=[dsiall.dsi_pval]';
osi_pval=[dsiall.osi_pval]';

% figure;
% tmp1=histcounts([dsiall(osi_pval<0.05).irat],.5:1:5.5);
% tmp2=histcounts([dsiall.irat],.5:1:5.5);
% subplot(121);
% histogram([dsiall(osi_pval<0.05).irat]);hold on;
% histogram([dsiall.irat]);
% xlabel('rat');ylabel('# of neuron');
% legend({'ori-selective','all'});
% subplot(122);
% tmpr=tmp1./tmp2;
% bar(tmpr*100);
% xlabel('rat');ylabel('prop of neuron(%)');

dsirand2=cat(2,dsiall.dsi_rand);
osirand2=cat(2,dsiall.osi_rand);
dsirandmean=mean(dsirand2,1)';
osirandmean=mean(osirand2,1)';
dsirand=dsirand2(:);
osirand=osirand2(:);

dsinorm=dsi./dsirandmean;
osinorm=osi./osirandmean;
dsi_fanofactorPval=cat(1,dsiall.fanofactor_dsi_pval);
osi_fanofactorPval=cat(1,dsiall.fanofactor_osi_pval);
pval_dosinorm=signrank(dsinorm,osinorm);
% figure;histogram(mean(dsirand2,1))
% figure;plot_scatter([dsinorm(:),osinorm(:)])

[dosi_m,dosi_se]=mean_se([dsi,osi,dsinorm,osinorm]);
[dosirand_m,dosirand_se]=mean_se([dsirand(:),osirand(:)]);

% figure;
% histogram(dsi(osi_pval>0.05),0:.1:1);
% hold on;histogram(dsi(osi_pval>0.05& dsi_pval<0.05),0:.1:1);
% xlabel('Direction selectivity index');ylabel('# of cells');
% legend({'nonsig OSI','nonsig OSI & sig DSI'});
% % find(sig(1,:)&~sig(2,:))

% preferred orientation
pref=zeros(length(dsiall),1);
for i=1:length(dsiall)
if mean(dsiall(i).peakrateMean([1,3]))>mean(dsiall(i).peakrateMean([2,4]))
    pref(i)=0; % vertical
else
    pref(i)=90;
end
end
histcounts(pref(osi_pval<0.05),[-1,1,89,91])

% figure;histogram(dsi_1track);
gooddsiind=find(~isnan(dsi_1track))';
f1=figure('position',[441   503   876   458]);
subplot(231);plot_scatter([dsi_1track(gooddsiind),osi(gooddsiind)],'.');xlabel('Bi-directionality');ylabel('OSI');
[r1,p1]=corr(dsi_1track(gooddsiind),osi(gooddsiind),'type','spearman');title(sprintf('r=%.2g,p=%.2g',r1,p1));
subplot(232);scatter(dsi_1track(gooddsiind),osinorm(gooddsiind),'.');xlabel('Bi-directionality');ylabel('OSI norm');
[r1,p1]=corr(dsi_1track(gooddsiind),osinorm(gooddsiind));title(sprintf('r=%.2g,p=%.2g',r1,p1));
subplot(233);scatter(dsi_1track(gooddsiind),osi_pval(gooddsiind),'.');hold on;plot([0,1],0.05*[1,1],'r');
[r1,p1]=corr(dsi_1track(gooddsiind),osi_pval(gooddsiind));title(sprintf('r=%.2g,p=%.2g',r1,p1));
xlabel('Bi-directionality');ylabel('OSI p-value');set(gca,'yscale','log');
subplot(234);plot_scatter([dsi_1track(gooddsiind),dsi(gooddsiind)],'.');xlabel('Bi-directionality');ylabel('DSI');
[r1,p1]=corr(dsi_1track(gooddsiind),dsi(gooddsiind));title(sprintf('r=%.2g,p=%.2g',r1,p1));
subplot(235);scatter(dsi_1track(gooddsiind),dsinorm(gooddsiind),'.');xlabel('Bi-directionality');ylabel('DSI norm');
[r1,p1]=corr(dsi_1track(gooddsiind),dsinorm(gooddsiind));title(sprintf('r=%.2g,p=%.2g',r1,p1));
subplot(236);scatter(dsi_1track(gooddsiind),dsi_pval(gooddsiind),'.');hold on;plot([0,1],0.05*[1,1],'r');
xlabel('Bi-directionality');ylabel('DSI p-value');set(gca,'yscale','log');
[r1,p1]=corr(dsi_1track(gooddsiind),dsi_pval(gooddsiind));title(sprintf('r=%.2g,p=%.2g',r1,p1));

setaxisformalAll


sig=[dsi_pval'<0.05; osi_pval'<0.05];
sig_str={'dsi all','dsi only','dsi&osi','osi only','osi all'};
sig_n=[sum(sig(1,:)),sum(sig(1,:)&~sig(2,:)),sum(sig(1,:)&sig(2,:)),sum(~sig(1,:)&sig(2,:)),sum(sig(2,:))];
f2=figure('position',[121   483   702   622]);
subplot(221);
bar(sig_n);
set(gca,'xtick',1:length(sig_n),'xticklabel',sig_str,'XTickLabelRotation',30,...
    'box','off');
ylabel('# of cells');
subplot(2,2,2);
[vh,vs]=venn([sig_n(1),sig_n(5)],sig_n(3));
text(vs.ZoneCentroid(:,1)-.3,vs.ZoneCentroid(:,2),num2strcell(sig_n([2,4,3])));
legend(vh,{'Direction','Orientation'});
axis off
subplot(2,2,3);
pie([sig_n(5)/length(osi),1-sig_n(5)/length(osi)],{'Ori sig','non-sig'});
title(sprintf('%d / %d',sig_n(5),length(osi)));

axis off

setaxisformalAll
%
f3=figure('position',[  63   239   655   594]);
subplot(221)
errorbarformal(1:4,[dosi_m(1),dosirand_m(1),dosi_m(2),dosirand_m(2)],...
    [dosi_se(1),dosirand_se(1),dosi_se(2),dosirand_se(2)],...
    {'DSI','DSIrand','OSI','OSIrand'});
set(gca,'xticklabelrotation',30)
subplot(222)
errorbarformal(1:2,[dosi_m(3:4)],[dosi_se(3:4)],{'DSInorm','OSInorm'});hold on;
sigstar([1,2],pval_dosinorm);
subplot(223)
plot_scatter([dsi(:),osi(:)],'.');
xlabel('DSI');ylabel('OSI');
title(sprintf('p=%.2g',signrank(dsi,osi)))
subplot(224)
plot_scatter([dsinorm(:),osinorm(:)],'.');
xlabel('DSI norm');ylabel('OSI norm');
title(sprintf('p=%.2g',pval_dosinorm))
setaxisformalAll

%% fig s2b s2d orientation selective cell OSI
edge=0:.1:1;
edgex=edge2x(edge);
n_dsi=histcounts(dsi,edge,'normalization','probability');
n_dsi_sig=histcounts(dsi(dsi_pval<0.05),edge); n_dsi_sig=n_dsi_sig/length(dsi);
n_dsirand=histcounts(dsirand,edge,'normalization','probability');
n_osi=histcounts(osi,edge,'normalization','probability');hold on;
n_osi_sig=histcounts(osi(osi_pval<0.05),edge);n_osi_sig=n_osi_sig/length(osi);
n_osirand=histcounts(osirand,edge,'normalization','probability');

[h,p_osi]=kstest2(osi,osirand)
[h,p_dsi]=kstest2(dsi,dsirand)


f4=figure('position',[65   585   839   347]);
subplot(121)
h3=bar(edgex,n_dsirand*100,1,'edgecolor','none','facealpha',0.5,'facecolor','r');hold on;
h1=bar(edgex,n_dsi*100,1,'edgecolor','none','facealpha',0.5,'facecolor','b');hold on;
h2=bar(edgex,n_dsi_sig*100,1,'edgecolor','none','facealpha',0.5,'facecolor',[0,0.5,0.5]);

plot(edgex,n_dsi*100,'b','linewidth',2);hold on;
plot(edgex,n_dsi_sig*100,':','color',[0,0.5,0.5],'linewidth',2);
plot(edgex,n_dsirand*100,'r','linewidth',2);
% set([h1;h2;h3],'linewidth',2)
xlabel('Direction selectivity');ylabel('%cell');
title(sprintf('%.1f%%significant(%d/%d)',sum(dsi_pval<0.05)/length(dsi)*100,...
    sum(dsi_pval<0.05),length(dsi)));
legend([h1;h2;h3],{'All pyr cell','significant cell','Shuffle'},'box','off');
box off
subplot(122)
h3=bar(edgex,n_osirand*100,1,'edgecolor','none','facealpha',0.5,'facecolor','r');hold on;
h1=bar(edgex,n_osi*100,1,'edgecolor','none','facealpha',0.5,'facecolor','b');hold on;
h2=bar(edgex,n_osi_sig*100,1,'edgecolor','none','facealpha',0.5,'facecolor',[0,0.5,0.5]);
plot(edgex,n_osi*100,'b','linewidth',2);hold on;
plot(edgex,n_osi_sig*100,':','color',[0,0.5,0.5],'linewidth',2);
plot(edgex,n_osirand*100,'r','linewidth',2);
% set([h1;h2;h3],'linewidth',2)
xlabel('Orientation selectivity');ylabel('%cell');
title(sprintf('%.1f%%significant(%d/%d)',sum(osi_pval<0.05)/length(osi)*100,...
    sum(osi_pval<0.05),length(osi)));
legend([h1;h2;h3],{'All pyr cell','significant cell','Shuffle'},'box','off');
box off
setaxisformalAll
savepath='/media/Big1/15tracks/figures/';

% savepdf(f1,fullfile(savepath,'biDirectionality_OSI.pdf'));
% savepdf(f2,fullfile(savepath,'Neurons_circleven_DSI_OSI.pdf'));
% savepdf(f3,fullfile(savepath,'compare_DSI_OSI.pdf'));
% savepdf(f4,fullfile(savepath,'OSI_DSI_selectivity_distribution.pdf'));

%% fig 2g s2c orientation selective cell 
plotind=[91,238,269,124,30];
addbase=10^-4;
f5=figure('position',[152          50         851         733]);
subplot(221)
scatter(osi+addbase,osi_pval+addbase)
hold on;
plot([addbase,1],0.05*[1,1],'r');
plot(osi(plotind)+addbase,osi_pval(plotind)+addbase,'*');
% plot([addbase,1],0.01*[1,1],'color',[1,0.5,0]);
set(gca,'xscale','linear','yscale','log','box','off','tickdir','out','ylim',...
    [10^-5,2],'xlim',[0,1]);
xlabel('Orientation selectivity');
ylabel('p value');

subplot(223)
scatter(osi+addbase,osi_fanofactorPval +addbase)
set(gca,'xscale','linear','yscale','linear','box','off','tickdir','out');
xlabel('Orientation selectivity');
ylabel('reliability(1-fanofactor/shuffle)');

[x,y]=remove_nan(osi,osi_fanofactorPval);
[rfano,pfano]=corr(x,y)

subplot(224)
scatter(osi_pval+addbase,osi_fanofactorPval +addbase)
hold on;
plot(0.05*[1,1],[addbase,1],'r');
set(gca,'xscale','log','yscale','linear','box','off','tickdir','out','xlim',...
    [10^-5,2]);
xlabel('p value(osi)');
ylabel('reliability(1-fanofactor/shuffle)');
setaxisformalAll

% savepdf(f5,fullfile(savepath,'OSI_and_reliability.pdf'));
%% fig s2g  s2i properties of ori select neuron  
clear
load(['/media/Big1/15tracks/HeadDirection/property_of_ori_neuron.mat'],'M1_note_param','M2_note_param','cellidall',...
    'fieldlengthThr','goodcell','goodcellinfo','goodind','info','m1','m1mean',...
    'm1mean_ori_compare','m1mean_ori_meanse','m1mean_ori_pval','m1plot','m2','m2mean',...
    'm2mean_ori_compare','m2mean_ori_meanse','m2mean_ori_pval','minrate','ori_pval',...
    'osi','prefOri','ratnames','singleTrSes','tagstr');

table([M1_note_param';M2_note_param'],[m1mean_ori_pval;m2mean_ori_pval])


m2plot=[2,3,4,6,7,8,9];m1plot=[1,2,3,4,5,6];
% m2plot=1:length(M2_note_param);m1plot=1:length(M1_note_param);
nplot=length(m1plot)+length(m2plot);
barcolor=[    0.5,0.3,1
    1,0.4,0.3
    ];


%  ------------- cdf plot
f1=figure('color','w','position',[ 205          67        1026         862]);
nrow=ceil(sqrt(nplot+1));
for i=1:length(m2plot)
    subplot(nrow,nrow,i);
    %    h1=cdfplot(m2mean_ori_compare{m2plot(i),1}); hold on; set(h1,'Color',barcolor(1,:));
    %    h2=cdfplot(m2mean_ori_compare{m2plot(i),2}); set(h2,'Color',barcolor(2,:));
    
    [y,x]=ecdf(m2mean_ori_compare{m2plot(i),1});
    h1=plot(x,y,'Color',barcolor(1,:)); hold on;
    m=m2mean_ori_meanse(m2plot(i),1,1); [~,tmpind]=min(abs(x-m));
    plot(m*[1,1],[0,y(tmpind(1))],'-','Color',barcolor(1,:));
    
    [y,x]=ecdf(m2mean_ori_compare{m2plot(i),2});
    h2=plot(x,y,'Color',barcolor(2,:));
    m=m2mean_ori_meanse(m2plot(i),2,1);   [~,tmpind]=min(abs(x-m));
    plot(m*[1,1],[0,y( tmpind(1) )],'-','Color',barcolor(2,:));
    
    xlabel(M2_note_param{m2plot(i)}); ylabel('prop of neurons');
    title(sprintf('p=%.2g',m2mean_ori_pval(m2plot(i))));
    set(gca,'xlim',[min(x) max(x)]);
    if strcmp(M2_note_param{m2plot(i)},'VarianceOfFieldRate')
        set(gca,'xlim',[0,50]);
    end
    %    set(gca,'XTickLabelRotation',30);
end
for i=1:length(m1plot)
    subplot(nrow,nrow,length(m2plot)+i);
    %    h1=cdfplot(m1mean_ori_compare{m1plot(i),1}); hold on; set(h1,'Color',barcolor(1,:));
    %    h2=cdfplot(m1mean_ori_compare{m1plot(i),2}); set(h2,'Color',barcolor(2,:));
    [y,x]=ecdf(m1mean_ori_compare{m1plot(i),1});
    h1=plot(x,y,'Color',barcolor(1,:)); hold on;
    m=m1mean_ori_meanse(m1plot(i),1,1);   [~,tmpind]=min(abs(x-m));
    plot(m*[1,1],[0,y( tmpind(1) )],'-','Color',barcolor(1,:));
    
    [y,x]=ecdf(m1mean_ori_compare{m1plot(i),2});
    h2=plot(x,y,'Color',barcolor(2,:));
    m=m1mean_ori_meanse(m1plot(i),2,1);   [~,tmpind]=min(abs(x-m));
    plot(m*[1,1],[0,y(tmpind(1) )],'-','Color',barcolor(2,:));
    set(gca,'xlim',[min(x) max(x)]);
    xlabel(M1_note_param{m1plot(i)}); ylabel('prop of neurons');
    title(sprintf('p=%.2g',m1mean_ori_pval(m1plot(i))));
end
subplot(nrow,nrow,nplot+1);
plot(rand(1))
axis off;
legend([h1,h2],{'Ori significant','non-sig'},'position',[0.76,0.13,0.14,0.05]);
setaxisformalAll(10)
% ca1 ca2 ca3 tetrode location in hippocampus
ori_neuron_ratio=zeros(5,1);
for irat=1:5
    ori_neuron_ratio(irat)=sum([info.irat]==irat &[info.oriSelective])/sum([info.irat]==irat);
end

order=[3,4,1,2,5];
f2=figure;
bar(ori_neuron_ratio(order));
set(gca,'xtick',1:5,'xticklabel',[]);xlabel('Anterior -> Posterior');
ylabel('proportion of ori selective neurons');
setaxisformal
f3=figure;
plot(1:5,ori_neuron_ratio(order),'o-');
set(gca,'xtick',1:5,'xticklabel',[]);xlabel('Tetrode location(Anterior -> Posterior)');
ylabel('proportion of ori selective neurons');xlim([0,6])
setaxisformal
savepath='/media/Big1/15tracks/figures/';
% savepdf(f1,fullfile(savepath,'ori_neuron_properties.pdf'));
% savepdf(f2,fullfile(savepath,'ori_neuron_tetrode_anterior2posterior.pdf'));
savepdf(f3,fullfile(savepath,'ori_neuron_tetrode_anterior2posterior_dot.pdf'));
 %% fig s2h; tetrode location vs prop of ori cell
 clear
load('/media/Big1/15tracks/TetrodeLocation/tetrodelocation.mat');

 f1=figure('position',[131         138        861         440],'color','w');
ax1=subplot(1,3,[2,3]);
xyc=[];
for irat=1:5
   thisrat=find([ttinfo.irat]==irat & [ttinfo.nUnit]>=3);
   xyc=[xyc;  abs([ttinfo(thisrat).xmm] + drivecentermm(irat,1))',[ttinfo(thisrat).ymm]' + drivecentermm(irat,2),[ttinfo(thisrat).propOfOriSig]'];
   
   h1=scatter(abs([ttinfo(thisrat).xmm] + drivecentermm(irat,1)),[ttinfo(thisrat).ymm] + drivecentermm(irat,2),...
      46,[ttinfo(thisrat).propOfOriSig],'filled');hold on;
   set(gca,'ydir','reverse','xtick',0:10,'ytick',0:.5:10);
%    axis tight equal;
   ylabel('Posterior <-> Anterior (mm)');xlabel('Mid <-> Lateral (mm)');
end
colorbar;
yedge=linspace(min(xyc(:,2)),max(xyc(:,2)),5); yedgex=edge2x(yedge)';
[~,~,ia]=histcounts(xyc(:,2),yedge);
[cym,cyse]=splitapply(@mean_se,xyc(:,3),ia);
ax2=subplot(1,3,1);
% errorbar(cym,yedgex,zeros(size(yedgex)),zeros(size(yedgex)),cyse,cyse,'linewidth',2)
errorshade21(cym,yedgex,zeros(size(yedgex)),cyse,{'b.-','linewidth',2},{'b','FaceAlpha',0.5})
linkaxes([ax1,ax2],'y');
set(ax2,'ylim',get(ax1,'ylim'),'ydir','reverse','xlim',[0,0.5],'xtick',0:0.25:1,'ytick',0:.5:10);
ylabel('Posterior <-> Anterior (mm)');xlabel('Prop. of Ori. neuron');
% ax2.Position=[ax2.Position(1),ax1.Position(2),ax2.Position(3),ax1.Position(4)];
% ax3=subplot(2,1,2);
% shadedErrorBar(yedgex,cym,cyse,{'b'});
% set(ax3,'xlim',get(ax1,'ylim'),'ylim',[0,0.5]);
load colormapCopperLight.mat;
set(gcf,'colormap',cmap);
setaxisformalAll(14,0)

savepath='/media/Big1/15tracks/figures/';
% savepdf(f1,fullfile(savepath,'tetrode_location_prop_of_ori_cell.pdf'));
